import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
import pickle
import arviz as az
import pymc3 as pm
WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.
from matplotlib.colors import to_rgb
import scipy.stats as stats
from IPython.display import display
import matplotlib as mpl
%load_ext autoreload
%autoreload 2
import plotting_lib
writeOut = True
outPathPlots = "../plots/statistical_model_two_factors/"
outPathData = "../derived_data/statistical_model_two_factors/"
widthMM = 190
widthInch = widthMM / 25.4
ratio = 0.66666
heigthInch = ratio*widthInch
SMALL_SIZE = 8
MEDIUM_SIZE = 10
BIGGER_SIZE = 12
plt.rc('font', size=SMALL_SIZE) # controls default text sizes
plt.rc('axes', titlesize=SMALL_SIZE) # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title
sns.set_style("ticks")
dpi = 300
sizes = [SMALL_SIZE,MEDIUM_SIZE,BIGGER_SIZE]
numSamples = 500
numCores = 10
numTune = 1000
numPredSamples = 2000
random_seed=3651435
datafile = "../derived_data/preprocessing/preprocessed.dat"
with open(datafile, "rb") as f:
x1,x2,_,df,dataZ,dictMeanStd,dictTreatment,dictSoftware = pickle.load(f)
Show that everything is correct:
display(pd.DataFrame.from_dict({'x1':x1,'x2':x2}))
| x1 | x2 | |
|---|---|---|
| 0 | 0 | 5 |
| 1 | 1 | 5 |
| 2 | 0 | 5 |
| 3 | 1 | 5 |
| 4 | 0 | 5 |
| ... | ... | ... |
| 275 | 1 | 9 |
| 276 | 0 | 9 |
| 277 | 1 | 9 |
| 278 | 0 | 9 |
| 279 | 1 | 9 |
280 rows × 2 columns
x1 indicates the software used, x2 indicates the treatment applied.
for surfaceParam,(mean,std) in dictMeanStd.items():
print("Surface parameter {} has mean {} and standard deviation {}".format(surfaceParam,mean,std))
Surface parameter epLsar has mean 0.002866777015225286 and standard deviation 0.0019173233323041528 Surface parameter R² has mean 0.9973197285182144 and standard deviation 0.006745323717352575 Surface parameter Asfc has mean 16.912804592866785 and standard deviation 16.042490777228107 Surface parameter Smfc has mean 2.589874101745179 and standard deviation 10.663178442785044 Surface parameter HAsfc9 has mean 0.31872043020221136 and standard deviation 0.22913790943445264 Surface parameter HAsfc81 has mean 0.5885203613280154 and standard deviation 0.42897734163543366
for k,v in dictTreatment.items():
print("Number {} encodes treatment {}".format(k,v))
Number 5 encodes treatment Dry bamboo Number 6 encodes treatment Dry grass Number 7 encodes treatment Dry lucerne Number 0 encodes treatment BrushDirt Number 1 encodes treatment BrushNoDirt Number 4 encodes treatment Control Number 10 encodes treatment RubDirt Number 2 encodes treatment Clover Number 3 encodes treatment Clover+dust Number 8 encodes treatment Grass Number 9 encodes treatment Grass+dust
for k,v in dictSoftware.items():
print("Number {} encodes software {}".format(k,v))
Number 0 encodes software ConfoMap Number 1 encodes software Toothfrax
display(dataZ)
| index | TreatmentNumber | SoftwareNumber | DatasetNumber | NameNumber | epLsar_z | R²_z | Asfc_z | Smfc_z | HAsfc9_z | HAsfc81_z | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 5 | 0 | 0 | 116 | 0.896441 | 0.202280 | -0.548753 | -0.173680 | -0.872185 | -0.512580 |
| 1 | 1 | 5 | 1 | 0 | 116 | 0.967089 | 0.332122 | -0.410913 | -0.231700 | -0.799735 | -0.528438 |
| 2 | 2 | 5 | 0 | 0 | 117 | 1.663582 | 0.205805 | -0.406276 | -0.184263 | -0.625739 | -0.636413 |
| 3 | 3 | 5 | 1 | 0 | 117 | 1.559060 | 0.318335 | -0.231484 | -0.231700 | -0.652395 | -0.762985 |
| 4 | 4 | 5 | 0 | 0 | 118 | 1.235447 | 0.140997 | -0.524577 | -0.187419 | -0.425552 | -0.460441 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 275 | 275 | 9 | 1 | 2 | 51 | 0.812186 | 0.302472 | -0.940814 | -0.110773 | 2.259947 | 1.219611 |
| 276 | 276 | 9 | 0 | 2 | 52 | 0.279516 | 0.134624 | -0.883003 | -0.199355 | 1.616592 | 2.457881 |
| 277 | 277 | 9 | 1 | 2 | 52 | 0.141981 | 0.358659 | -0.882314 | -0.230373 | 2.779893 | 2.898056 |
| 278 | 278 | 9 | 0 | 2 | 53 | -0.859089 | 0.269132 | -0.964638 | -0.118959 | 0.316720 | 0.540611 |
| 279 | 279 | 9 | 1 | 2 | 53 | -1.128540 | 0.376153 | -0.964978 | -0.204577 | 0.119539 | 0.586105 |
280 rows × 11 columns
display(df)
| Dataset | Name | Software | Diet | Treatment | Before.after | epLsar | R² | Asfc | Smfc | HAsfc9 | HAsfc81 | NewEplsar | TreatmentNumber | SoftwareNumber | DatasetNumber | NameNumber | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | GuineaPigs | capor_2CC6B1_txP4_#1_1_100xL_1 | ConfoMap | Dry bamboo | Dry bamboo | NaN | 0.004586 | 0.998684 | 8.109445 | 0.737898 | 0.118870 | 0.368635 | 0.019529 | 5 | 0 | 0 | 116 |
| 1 | GuineaPigs | capor_2CC6B1_txP4_#1_1_100xL_1 | Toothfrax | Dry bamboo | Dry bamboo | NaN | 0.004721 | 0.999560 | 10.320730 | 0.119219 | 0.135471 | 0.361833 | NaN | 5 | 1 | 0 | 116 |
| 2 | GuineaPigs | capor_2CC6B1_txP4_#1_1_100xL_2 | ConfoMap | Dry bamboo | Dry bamboo | NaN | 0.006056 | 0.998708 | 10.395128 | 0.625040 | 0.175340 | 0.315513 | 0.020162 | 5 | 0 | 0 | 117 |
| 3 | GuineaPigs | capor_2CC6B1_txP4_#1_1_100xL_2 | Toothfrax | Dry bamboo | Dry bamboo | NaN | 0.005856 | 0.999467 | 13.199232 | 0.119219 | 0.169232 | 0.261217 | NaN | 5 | 1 | 0 | 117 |
| 4 | GuineaPigs | capor_2CC6B1_txP4_#1_1_100xL_3 | ConfoMap | Dry bamboo | Dry bamboo | NaN | 0.005236 | 0.998271 | 8.497286 | 0.591396 | 0.221210 | 0.391002 | 0.019804 | 5 | 0 | 0 | 118 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 275 | Sheeps | L8-Ovis-90730-lm2sin-a | Toothfrax | Grass+dust | Grass+dust | NaN | 0.004424 | 0.999360 | 1.819802 | 1.408678 | 0.836560 | 1.111706 | NaN | 9 | 1 | 2 | 51 |
| 276 | Sheeps | L8-Ovis-90764-lm2sin-a | ConfoMap | Grass+dust | Grass+dust | NaN | 0.003403 | 0.998228 | 2.747241 | 0.464115 | 0.689143 | 1.642896 | 0.018978 | 9 | 0 | 2 | 52 |
| 277 | Sheeps | L8-Ovis-90764-lm2sin-a | Toothfrax | Grass+dust | Grass+dust | NaN | 0.003139 | 0.999739 | 2.758297 | 0.133366 | 0.955699 | 1.831721 | NaN | 9 | 1 | 2 | 52 |
| 278 | Sheeps | L8-Ovis-90814-lm2sin-a | ConfoMap | Grass+dust | Grass+dust | NaN | 0.001220 | 0.999135 | 1.437607 | 1.321398 | 0.391293 | 0.820430 | 0.017498 | 9 | 0 | 2 | 53 |
| 279 | Sheeps | L8-Ovis-90814-lm2sin-a | Toothfrax | Grass+dust | Grass+dust | NaN | 0.000703 | 0.999857 | 1.432148 | 0.408433 | 0.346111 | 0.839946 | NaN | 9 | 1 | 2 | 53 |
280 rows × 17 columns
class TwoFactorModel(pm.Model):
"""
Compute params of priors and hyperpriors.
"""
def getParams(self,x1,x2,y):
# get lengths
Nx1Lvl = np.unique(x1).size
Nx2Lvl = np.unique(x2).size
dims = (Nx1Lvl, Nx2Lvl)
### get standard deviations
# convert to pandas dataframe to use their logic
df = pd.DataFrame.from_dict({'x1':x1,'x2':x2,'y':y})
s1 = df.groupby('x1').std()['y'].max()
s2 = df.groupby('x2').std()['y'].max()
stdSingle = (s1, s2)
prefac = 0.05
s12 = prefac * np.linalg.norm([s1,s2])
stdMulti = (s12)
return (dims, stdSingle, stdMulti)
def printParams(self,x1,x2,y):
dims, stdSingle, stdMulti = self.getParams(x1,x2,y)
Nx1Lvl, Nx2Lvl = dims
s1, s2 = stdSingle
s12 = stdMulti
print("The number of levels of the x variables are {}".format(dims))
print("The standard deviations used for the beta priors are {}".format(stdSingle))
print("The standard deviations used for the M12 priors are {}".format(stdMulti))
def __init__(self,name,x1,x2,y,model=None):
# call super's init first, passing model and name
super().__init__(name, model)
# get parameter of hyperpriors
dims, stdSingle, stdMulti = self.getParams(x1,x2,y)
Nx1Lvl, Nx2Lvl = dims
s1, s2 = stdSingle
s12 = stdMulti
### hyperpriors ###
# observation hyperpriors
lamY = 1/30.
muGamma = 0.5
sigmaGamma = 2.
# prediction hyperpriors
sigma0 = pm.HalfNormal('sigma0',sd=1)
sigma1 = pm.HalfNormal('sigma1',sd=s1, shape=Nx1Lvl)
sigma2 = pm.HalfNormal('sigma2',sd=s2, shape=Nx2Lvl)
beta2 = (np.sqrt(6)*sigma2)/(np.pi)
mu_b0 = pm.Normal('mu_b0', mu=0., sd=1)
mu_b1 = pm.Normal('mu_b1', mu=0., sd=1, shape=Nx1Lvl)
mu_b2 = pm.Normal('mu_b2', mu=0., sd=1, shape=Nx2Lvl)
sigma12 = pm.HalfNormal('sigma12',sd=s12)
### priors ###
# observation priors
nuY = pm.Exponential('nuY',lam=lamY)
sigmaY = pm.Gamma('sigmaY',mu=muGamma, sigma=sigmaGamma)
# prediction priors
b0_dist = pm.Normal('b0_dist', mu=0, sd=1)
b0 = pm.Deterministic("b0", mu_b0 + b0_dist * sigma0)
b1_dist = pm.Normal('b1_dist', mu=0, sd=1)
b1 = pm.Deterministic("b1", mu_b1 + b1_dist * sigma1)
b2_beta = pm.HalfNormal('b2_beta', sd=beta2, shape=Nx2Lvl)
b2_dist = pm.Gumbel('b2_dist', mu=0, beta=1)
b2 = pm.Deterministic("b2", mu_b2 + b2_beta * b2_dist)
mu_M12 = pm.Normal('mu_M12', mu=0., sd=1, shape=[Nx1Lvl, Nx2Lvl])
M12_dist = pm.Normal('M12_dist', mu=0, sd=1)
M12 = pm.Deterministic("M12", mu_M12 + M12_dist * sigma12)
#### prediction ###
mu = pm.Deterministic('mu',b0 + b1[x1]+ b2[x2] + M12[x1,x2] )
### observation ###
y = pm.StudentT('y',nu = nuY, mu=mu, sd=sigmaY, observed=y)
with pm.Model() as model:
epLsarModel = TwoFactorModel('epLsar',x1,x2,dataZ.epLsar_z.values)
epLsarModel.printParams(x1,x2,dataZ.epLsar_z.values)
The number of levels of the x variables are (2, 11) The standard deviations used for the beta priors are (1.0100483277420549, 1.278487256789849) The standard deviations used for the M12 priors are 0.08146666941376325
try:
graph_epLsar = pm.model_to_graphviz(epLsarModel)
except:
graph_epLsar = "Could not make graph"
graph_epLsar
with epLsarModel as model:
prior_pred_epLsar = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)
plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_epLsar,dataZ.epLsar_z.values,'epLsar')
Prior choice is as intended: Broad over the data range.
with epLsarModel as model:
trace_epLsar = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.95,random_seed=random_seed)
#fit_epLsar = pm.fit(random_seed=random_seed)
#trace_epLsar = fit_epLsar.sample(draws=numSamples)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (10 chains in 10 jobs) NUTS: [epLsar_M12_dist, epLsar_mu_M12, epLsar_b2_dist, epLsar_b2_beta, epLsar_b1_dist, epLsar_b0_dist, epLsar_sigmaY, epLsar_nuY, epLsar_sigma12, epLsar_mu_b2, epLsar_mu_b1, epLsar_mu_b0, epLsar_sigma2, epLsar_sigma1, epLsar_sigma0]
Sampling 10 chains for 1_000 tune and 500 draw iterations (10_000 + 5_000 draws total) took 183 seconds.
with epLsarModel as model:
if writeOut:
with open(outPathData + 'model_{}.pkl'.format('epLsar'), 'wb') as buff:
pickle.dump({'model':epLsarModel, 'trace': trace_epLsar}, buff)
if writeOut:
np.save('../derived_data/statistical_model_two_factors/epLsar_oldb1', trace_epLsar['epLsar_b1'])
np.save('../derived_data/statistical_model_two_factors/epLsar_oldb2', trace_epLsar['epLsar_b2'])
np.save('../derived_data/statistical_model_two_factors/epLsar_oldM12', trace_epLsar['epLsar_M12'])
with epLsarModel as model:
dataTrace_epLsar = az.from_pymc3(trace=trace_epLsar)
pm.summary(dataTrace_epLsar,hdi_prob=0.95).round(2)
| mean | sd | hdi_2.5% | hdi_97.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| epLsar_mu_b0 | -0.11 | 0.81 | -1.67 | 1.48 | 0.01 | 0.01 | 5029.0 | 2319.0 | 5029.0 | 3627.0 | 1.0 |
| epLsar_mu_b1[0] | 0.00 | 0.73 | -1.48 | 1.36 | 0.01 | 0.01 | 4766.0 | 3057.0 | 4763.0 | 3988.0 | 1.0 |
| epLsar_mu_b1[1] | -0.07 | 0.71 | -1.47 | 1.34 | 0.01 | 0.01 | 4575.0 | 2861.0 | 4587.0 | 3606.0 | 1.0 |
| epLsar_mu_b2[0] | -0.25 | 0.66 | -1.63 | 0.94 | 0.01 | 0.01 | 3705.0 | 2967.0 | 3733.0 | 4054.0 | 1.0 |
| epLsar_mu_b2[1] | -0.36 | 0.67 | -1.64 | 0.98 | 0.01 | 0.01 | 4093.0 | 3216.0 | 4098.0 | 3736.0 | 1.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| epLsar_mu[275] | 0.66 | 0.27 | 0.12 | 1.16 | 0.00 | 0.00 | 5124.0 | 4919.0 | 5132.0 | 4303.0 | 1.0 |
| epLsar_mu[276] | 0.68 | 0.31 | 0.09 | 1.30 | 0.00 | 0.00 | 5093.0 | 5024.0 | 5103.0 | 3850.0 | 1.0 |
| epLsar_mu[277] | 0.66 | 0.27 | 0.12 | 1.16 | 0.00 | 0.00 | 5124.0 | 4919.0 | 5132.0 | 4303.0 | 1.0 |
| epLsar_mu[278] | 0.68 | 0.31 | 0.09 | 1.30 | 0.00 | 0.00 | 5093.0 | 5024.0 | 5103.0 | 3850.0 | 1.0 |
| epLsar_mu[279] | 0.66 | 0.27 | 0.12 | 1.16 | 0.00 | 0.00 | 5124.0 | 4919.0 | 5132.0 | 4303.0 | 1.0 |
384 rows × 11 columns
plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_epLsar,dataTrace_epLsar,'epLsar')
/home/bob/.local/lib/python3.8/site-packages/arviz/plots/backends/matplotlib/pairplot.py:212: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid warnings.warn(
with epLsarModel as model:
plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_epLsar,'epLsar')
with epLsarModel as model:
plotting_lib.pm.energyplot(trace_epLsar)
with epLsarModel as model:
posterior_pred_epLsar = pm.sample_posterior_predictive(trace_epLsar,samples=numPredSamples,random_seed=random_seed)
/home/bob/.local/lib/python3.8/site-packages/pymc3/sampling.py:1707: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample warnings.warn(
plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_epLsar,posterior_pred_epLsar,dataZ.epLsar_z.values,'epLsar')
plotting_lib.plotLevels(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_epLsar,'epLsar',x1,x2)
plotting_lib.plotLevelsStd(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_epLsar,'epLsar',x1,x2)
df_hdi_epLsar = plotting_lib.plotTreatmentPosterior(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_epLsar,'epLsar',x1,x2)
df_hdi_epLsar
| Treatment_i | Treatment_j | hdi_ConfoMap_2.5% | hdi_ConfoMap_97.5% | isSignificant_on_ConfoMap | hdi_Toothfrax_2.5% | hdi_Toothfrax_97.5% | isSignificant_on_Toothfrax | |
|---|---|---|---|---|---|---|---|---|
| 0 | Dry grass | Dry bamboo | -0.004200 | -0.002884 | True | -0.003967 | -0.002636 | True |
| 1 | Dry lucerne | Dry bamboo | -0.003305 | -0.001994 | True | -0.003411 | -0.002161 | True |
| 2 | Dry lucerne | Dry grass | 0.000228 | 0.001524 | True | -0.000167 | 0.001116 | False |
| 3 | BrushNoDirt | BrushDirt | -0.001227 | 0.001107 | False | -0.001528 | 0.000527 | False |
| 4 | Control | BrushDirt | -0.002046 | 0.000359 | False | -0.001954 | 0.000235 | False |
| 5 | Control | BrushNoDirt | -0.001872 | 0.000384 | False | -0.001414 | 0.000676 | False |
| 6 | RubDirt | BrushDirt | -0.001636 | 0.000627 | False | -0.001497 | 0.000683 | False |
| 7 | RubDirt | BrushNoDirt | -0.001442 | 0.000661 | False | -0.000890 | 0.001130 | False |
| 8 | RubDirt | Control | -0.000753 | 0.001508 | False | -0.000598 | 0.001547 | False |
| 9 | Clover+dust | Clover | -0.001140 | 0.001047 | False | -0.000746 | 0.001576 | False |
| 10 | Grass | Clover | -0.000153 | 0.002973 | False | 0.000293 | 0.003795 | True |
| 11 | Grass | Clover+dust | -0.000151 | 0.002972 | False | -0.000175 | 0.003389 | False |
| 12 | Grass+dust | Clover | 0.000533 | 0.003296 | True | 0.000819 | 0.003383 | True |
| 13 | Grass+dust | Clover+dust | 0.000445 | 0.003307 | True | 0.000379 | 0.002978 | True |
| 14 | Grass+dust | Grass | -0.001303 | 0.002230 | False | -0.001721 | 0.001958 | False |
if writeOut:
df_hdi_epLsar.to_csv(outPathData+ 'hdi_{}.csv'.format('epLsar'))
with pm.Model() as model:
RsquaredModel = TwoFactorModel('R²',x1,x2,dataZ["R²_z"].values)
RsquaredModel.printParams(x1,x2,dataZ["R²_z"].values)
The number of levels of the x variables are (2, 11) The standard deviations used for the beta priors are (1.3433087377122408, 2.3171005331515255) The standard deviations used for the M12 priors are 0.13391632877981252
pm.model_to_graphviz(RsquaredModel)
with RsquaredModel as model:
prior_pred_Rsquared = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)
plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Rsquared,dataZ["R²_z"].values,'R²')
with RsquaredModel as model:
trace_Rsquared = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.99,random_seed=random_seed)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (10 chains in 10 jobs) NUTS: [R²_M12_dist, R²_mu_M12, R²_b2_dist, R²_b2_beta, R²_b1_dist, R²_b0_dist, R²_sigmaY, R²_nuY, R²_sigma12, R²_mu_b2, R²_mu_b1, R²_mu_b0, R²_sigma2, R²_sigma1, R²_sigma0]
Sampling 10 chains for 1_000 tune and 500 draw iterations (10_000 + 5_000 draws total) took 3971 seconds.
with RsquaredModel as model:
if writeOut:
with open(outPathData + 'model_{}.pkl'.format('R²'), 'wb') as buff:
pickle.dump({'model': RsquaredModel, 'trace': trace_Rsquared}, buff)
with RsquaredModel as model:
dataTrace_Rsquared = az.from_pymc3(trace=trace_Rsquared)
pm.summary(dataTrace_Rsquared,hdi_prob=0.95).round(2)
| mean | sd | hdi_2.5% | hdi_97.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| R²_mu_b0 | 0.07 | 0.83 | -1.55 | 1.64 | 0.01 | 0.01 | 7950.0 | 2327.0 | 7929.0 | 4058.0 | 1.0 |
| R²_mu_b1[0] | -0.05 | 0.73 | -1.49 | 1.35 | 0.01 | 0.01 | 6903.0 | 2379.0 | 6902.0 | 3964.0 | 1.0 |
| R²_mu_b1[1] | 0.11 | 0.72 | -1.23 | 1.58 | 0.01 | 0.01 | 7161.0 | 2924.0 | 7164.0 | 4300.0 | 1.0 |
| R²_mu_b2[0] | -0.01 | 0.66 | -1.24 | 1.29 | 0.01 | 0.01 | 5941.0 | 2544.0 | 5938.0 | 3984.0 | 1.0 |
| R²_mu_b2[1] | -0.05 | 0.66 | -1.31 | 1.29 | 0.01 | 0.01 | 6122.0 | 2330.0 | 6086.0 | 3681.0 | 1.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| R²_mu[275] | 0.34 | 0.01 | 0.31 | 0.36 | 0.00 | 0.00 | 4973.0 | 4968.0 | 4975.0 | 4789.0 | 1.0 |
| R²_mu[276] | 0.11 | 0.10 | -0.03 | 0.26 | 0.00 | 0.00 | 4537.0 | 4387.0 | 4704.0 | 4535.0 | 1.0 |
| R²_mu[277] | 0.34 | 0.01 | 0.31 | 0.36 | 0.00 | 0.00 | 4973.0 | 4968.0 | 4975.0 | 4789.0 | 1.0 |
| R²_mu[278] | 0.11 | 0.10 | -0.03 | 0.26 | 0.00 | 0.00 | 4537.0 | 4387.0 | 4704.0 | 4535.0 | 1.0 |
| R²_mu[279] | 0.34 | 0.01 | 0.31 | 0.36 | 0.00 | 0.00 | 4973.0 | 4968.0 | 4975.0 | 4789.0 | 1.0 |
384 rows × 11 columns
plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_Rsquared,dataTrace_Rsquared,'R²')
/home/bob/.local/lib/python3.8/site-packages/arviz/plots/backends/matplotlib/pairplot.py:212: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid warnings.warn(
with RsquaredModel as model:
plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_Rsquared,'R²')
with RsquaredModel as model:
plotting_lib.pm.energyplot(trace_Rsquared)
with RsquaredModel as model:
posterior_pred_Rsquared = pm.sample_posterior_predictive(trace_Rsquared,samples=numPredSamples,random_seed=random_seed)
/home/bob/.local/lib/python3.8/site-packages/pymc3/sampling.py:1707: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample warnings.warn(
plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Rsquared,posterior_pred_Rsquared,dataZ["R²_z"].values,'R²')
with RsquaredModel as model:
pm_data_Rsquared = az.from_pymc3(trace=trace_Rsquared,prior=prior_pred_Rsquared,posterior_predictive=posterior_pred_Rsquared)
arviz.data.io_pymc3 - WARNING - posterior predictive variable R²_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.
plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_Rsquared,'R²')
plotting_lib.plotLevels(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Rsquared,'R²',x1,x2)
plotting_lib.plotLevelsStd(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Rsquared,'R²',x1,x2)
plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_Rsquared,'R²')
df_hdi_R = plotting_lib.plotTreatmentPosterior(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Rsquared,'R²',x1,x2)
/home/bob/.local/lib/python3.8/site-packages/seaborn/distributions.py:305: UserWarning: Dataset has 0 variance; skipping density estimate. warnings.warn(msg, UserWarning)
df_hdi_R
| Treatment_i | Treatment_j | hdi_ConfoMap_2.5% | hdi_ConfoMap_97.5% | isSignificant_on_ConfoMap | hdi_Toothfrax_2.5% | hdi_Toothfrax_97.5% | isSignificant_on_Toothfrax | |
|---|---|---|---|---|---|---|---|---|
| 0 | Dry grass | Dry bamboo | -0.000285 | 0.000175 | False | -0.000257 | 0.000188 | False |
| 1 | Dry lucerne | Dry bamboo | -0.000116 | 0.000248 | False | -0.000326 | 0.000121 | False |
| 2 | Dry lucerne | Dry grass | -0.000115 | 0.000330 | False | -0.000349 | 0.000179 | False |
| 3 | BrushNoDirt | BrushDirt | -0.002096 | -0.000672 | True | -0.000420 | 0.000635 | False |
| 4 | Control | BrushDirt | -0.001249 | 0.000041 | False | -0.000416 | 0.000408 | False |
| 5 | Control | BrushNoDirt | -0.000080 | 0.001656 | False | -0.000801 | 0.000387 | False |
| 6 | RubDirt | BrushDirt | -0.001154 | 0.000799 | False | -0.000185 | 0.000635 | False |
| 7 | RubDirt | BrushNoDirt | 0.000181 | 0.002386 | True | -0.000556 | 0.000619 | False |
| 8 | RubDirt | Control | -0.000680 | 0.001489 | False | -0.000252 | 0.000722 | False |
| 9 | Clover+dust | Clover | -0.001005 | -0.000290 | True | -0.000138 | 0.000401 | False |
| 10 | Grass | Clover | -0.001250 | -0.000102 | True | -0.000719 | -0.000045 | True |
| 11 | Grass | Clover+dust | -0.000639 | 0.000580 | False | -0.000792 | -0.000199 | True |
| 12 | Grass+dust | Clover | -0.001970 | 0.000084 | False | -0.000435 | 0.000140 | False |
| 13 | Grass+dust | Clover+dust | -0.001314 | 0.000824 | False | -0.000515 | -0.000045 | True |
| 14 | Grass+dust | Grass | -0.001536 | 0.000874 | False | -0.000091 | 0.000540 | False |
if writeOut:
df_hdi_R.to_csv(outPathData+ 'hdi_{}.csv'.format('R²'))
with pm.Model() as model:
AsfcModel = TwoFactorModel('Asfc',x1,x2,dataZ["Asfc_z"].values)
AsfcModel.printParams(x1,x2,dataZ["Asfc_z"].values)
The number of levels of the x variables are (2, 11) The standard deviations used for the beta priors are (1.0620939450667206, 0.7885965314025212) The standard deviations used for the M12 priors are 0.06614242279897746
pm.model_to_graphviz(AsfcModel)
with AsfcModel as model:
prior_pred_Asfc = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)
plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Asfc,dataZ["Asfc_z"].values,'Asfc')
Prior choice is as intended: Broad over the data range.
with AsfcModel as model:
trace_Asfc = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.99,random_seed=random_seed)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (10 chains in 10 jobs) NUTS: [Asfc_M12_dist, Asfc_mu_M12, Asfc_b2_dist, Asfc_b2_beta, Asfc_b1_dist, Asfc_b0_dist, Asfc_sigmaY, Asfc_nuY, Asfc_sigma12, Asfc_mu_b2, Asfc_mu_b1, Asfc_mu_b0, Asfc_sigma2, Asfc_sigma1, Asfc_sigma0]
Sampling 10 chains for 1_000 tune and 500 draw iterations (10_000 + 5_000 draws total) took 1032 seconds.
with AsfcModel as model:
if writeOut:
with open(outPathData + 'model_{}.pkl'.format('Asfc'), 'wb') as buff:
pickle.dump({'model': AsfcModel, 'trace': trace_Asfc}, buff)
with AsfcModel as model:
dataTrace_Asfc = az.from_pymc3(trace=trace_Asfc)
pm.summary(dataTrace_Asfc,hdi_prob=0.95).round(2)
| mean | sd | hdi_2.5% | hdi_97.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Asfc_mu_b0 | 0.03 | 0.79 | -1.54 | 1.57 | 0.01 | 0.01 | 7171.0 | 2541.0 | 7170.0 | 3979.0 | 1.0 |
| Asfc_mu_b1[0] | -0.08 | 0.70 | -1.44 | 1.28 | 0.01 | 0.01 | 6377.0 | 2385.0 | 6384.0 | 3933.0 | 1.0 |
| Asfc_mu_b1[1] | 0.11 | 0.71 | -1.33 | 1.47 | 0.01 | 0.01 | 6789.0 | 2647.0 | 6777.0 | 4089.0 | 1.0 |
| Asfc_mu_b2[0] | 0.67 | 0.68 | -0.75 | 1.92 | 0.01 | 0.01 | 6024.0 | 4059.0 | 6090.0 | 3494.0 | 1.0 |
| Asfc_mu_b2[1] | 0.80 | 0.68 | -0.50 | 2.16 | 0.01 | 0.01 | 6900.0 | 4532.0 | 6872.0 | 3974.0 | 1.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| Asfc_mu[275] | -0.94 | 0.05 | -1.04 | -0.84 | 0.00 | 0.00 | 4903.0 | 4900.0 | 4919.0 | 4492.0 | 1.0 |
| Asfc_mu[276] | -0.94 | 0.05 | -1.04 | -0.84 | 0.00 | 0.00 | 5312.0 | 5312.0 | 5324.0 | 4413.0 | 1.0 |
| Asfc_mu[277] | -0.94 | 0.05 | -1.04 | -0.84 | 0.00 | 0.00 | 4903.0 | 4900.0 | 4919.0 | 4492.0 | 1.0 |
| Asfc_mu[278] | -0.94 | 0.05 | -1.04 | -0.84 | 0.00 | 0.00 | 5312.0 | 5312.0 | 5324.0 | 4413.0 | 1.0 |
| Asfc_mu[279] | -0.94 | 0.05 | -1.04 | -0.84 | 0.00 | 0.00 | 4903.0 | 4900.0 | 4919.0 | 4492.0 | 1.0 |
384 rows × 11 columns
plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_Asfc,dataTrace_Asfc,'Asfc')
/home/bob/.local/lib/python3.8/site-packages/arviz/plots/backends/matplotlib/pairplot.py:212: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid warnings.warn(
with AsfcModel as model:
plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_Asfc,'Asfc')
with AsfcModel as model:
plotting_lib.pm.energyplot(trace_Asfc)
with AsfcModel as model:
posterior_pred_Asfc = pm.sample_posterior_predictive(trace_Asfc,samples=numPredSamples,random_seed=random_seed)
/home/bob/.local/lib/python3.8/site-packages/pymc3/sampling.py:1707: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample warnings.warn(
plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Asfc,posterior_pred_Asfc,dataZ["Asfc_z"].values,'Asfc')
plotting_lib.plotLevels(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Asfc,'Asfc',x1,x2)
plotting_lib.plotLevelsStd(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Asfc,'Asfc',x1,x2)
with AsfcModel as model:
pm_data_Asfc = az.from_pymc3(trace=trace_Asfc,prior=prior_pred_Asfc,posterior_predictive=posterior_pred_Asfc)
arviz.data.io_pymc3 - WARNING - posterior predictive variable Asfc_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.
plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_Asfc,'Asfc')
plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_Asfc,'Asfc')
df_hdi_Asfc = plotting_lib.plotTreatmentPosterior(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Asfc,'Asfc',x1,x2)
df_hdi_Asfc
| Treatment_i | Treatment_j | hdi_ConfoMap_2.5% | hdi_ConfoMap_97.5% | isSignificant_on_ConfoMap | hdi_Toothfrax_2.5% | hdi_Toothfrax_97.5% | isSignificant_on_Toothfrax | |
|---|---|---|---|---|---|---|---|---|
| 0 | Dry grass | Dry bamboo | -3.939803 | 0.935695 | False | -4.447570 | 0.965258 | False |
| 1 | Dry lucerne | Dry bamboo | -4.031412 | 0.625569 | False | -5.003380 | 0.239099 | False |
| 2 | Dry lucerne | Dry grass | -2.410119 | 1.424301 | False | -2.515884 | 1.276464 | False |
| 3 | BrushNoDirt | BrushDirt | 0.013059 | 7.915624 | True | -2.631889 | 6.887819 | False |
| 4 | Control | BrushDirt | 3.910899 | 18.434847 | True | 1.604132 | 14.612596 | True |
| 5 | Control | BrushNoDirt | 0.022112 | 14.415237 | True | -0.616464 | 13.182578 | False |
| 6 | RubDirt | BrushDirt | -0.259987 | 20.831039 | False | 1.502426 | 13.522475 | True |
| 7 | RubDirt | BrushNoDirt | -4.336062 | 16.137641 | False | -1.231718 | 11.842114 | False |
| 8 | RubDirt | Control | -14.247123 | 11.427003 | False | -8.595219 | 7.100929 | False |
| 9 | Clover+dust | Clover | -4.902063 | 2.035175 | False | -4.562186 | 1.844004 | False |
| 10 | Grass | Clover | -5.948810 | 1.110713 | False | -5.876442 | 0.688654 | False |
| 11 | Grass | Clover+dust | -3.539383 | 1.255674 | False | -3.501566 | 1.339794 | False |
| 12 | Grass+dust | Clover | -6.781348 | 0.056399 | False | -6.135492 | 0.254948 | False |
| 13 | Grass+dust | Clover+dust | -4.430594 | 0.335530 | False | -3.919102 | 0.673998 | False |
| 14 | Grass+dust | Grass | -3.353905 | 1.617864 | False | -3.044589 | 1.686165 | False |
if writeOut:
df_hdi_Asfc.to_csv(outPathData+ 'hdi_{}.csv'.format('Asfc'))
with pm.Model() as model:
SmfcModel = TwoFactorModel('Smfc',x1,x2,dataZ.Smfc_z.values)
SmfcModel.printParams(x1,x2,dataZ.Smfc_z.values)
The number of levels of the x variables are (2, 11) The standard deviations used for the beta priors are (1.190539275484547, 2.890262579386322) The standard deviations used for the M12 priors are 0.15629300643560595
pm.model_to_graphviz(SmfcModel)
with SmfcModel as model:
prior_pred_Smfc = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)
plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Smfc,dataZ.Smfc_z.values,'Smfc')
Prior choice is as intended: Broad over the data range.
with SmfcModel as model:
trace_Smfc = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.99,random_seed=random_seed)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (10 chains in 10 jobs) NUTS: [Smfc_M12_dist, Smfc_mu_M12, Smfc_b2_dist, Smfc_b2_beta, Smfc_b1_dist, Smfc_b0_dist, Smfc_sigmaY, Smfc_nuY, Smfc_sigma12, Smfc_mu_b2, Smfc_mu_b1, Smfc_mu_b0, Smfc_sigma2, Smfc_sigma1, Smfc_sigma0]
--------------------------------------------------------------------------- KeyboardInterrupt Traceback (most recent call last) ~/.local/lib/python3.8/site-packages/pymc3/sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, callback, discard_tuned_samples, mp_ctx, pickle_backend, **kwargs) 1485 with sampler: -> 1486 for draw in sampler: 1487 trace = traces[draw.chain - chain] ~/.local/lib/python3.8/site-packages/pymc3/parallel_sampling.py in __iter__(self) 491 while self._active: --> 492 draw = ProcessAdapter.recv_draw(self._active) 493 proc, is_last, draw, tuning, stats, warns = draw ~/.local/lib/python3.8/site-packages/pymc3/parallel_sampling.py in recv_draw(processes, timeout) 351 pipes = [proc._msg_pipe for proc in processes] --> 352 ready = multiprocessing.connection.wait(pipes) 353 if not ready: /usr/lib/python3.8/multiprocessing/connection.py in wait(object_list, timeout) 930 while True: --> 931 ready = selector.select(timeout) 932 if ready: /usr/lib/python3.8/selectors.py in select(self, timeout) 414 try: --> 415 fd_event_list = self._selector.poll(timeout) 416 except InterruptedError: KeyboardInterrupt: During handling of the above exception, another exception occurred: ValueError Traceback (most recent call last) <ipython-input-90-4de91012caf3> in <module> 1 with SmfcModel as model: ----> 2 trace_Smfc = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.99,random_seed=random_seed) ~/.local/lib/python3.8/site-packages/pymc3/sampling.py in sample(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, callback, return_inferencedata, idata_kwargs, mp_ctx, pickle_backend, **kwargs) 543 _print_step_hierarchy(step) 544 try: --> 545 trace = _mp_sample(**sample_args, **parallel_args) 546 except pickle.PickleError: 547 _log.warning("Could not pickle model, sampling singlethreaded.") ~/.local/lib/python3.8/site-packages/pymc3/sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, callback, discard_tuned_samples, mp_ctx, pickle_backend, **kwargs) 1510 except KeyboardInterrupt: 1511 if discard_tuned_samples: -> 1512 traces, length = _choose_chains(traces, tune) 1513 else: 1514 traces, length = _choose_chains(traces, 0) ~/.local/lib/python3.8/site-packages/pymc3/sampling.py in _choose_chains(traces, tune) 1528 lengths = [max(0, len(trace) - tune) for trace in traces] 1529 if not sum(lengths): -> 1530 raise ValueError("Not enough samples to build a trace.") 1531 1532 idxs = np.argsort(lengths)[::-1] ValueError: Not enough samples to build a trace.
Analysis stopped here because sampling did not converge. As the plot shows, some data points are very far away from the others, which would require the analysis to be based on more heavy-tailed distributions.
with pm.Model() as model:
HAsfc9Model = TwoFactorModel('HAsfc9',x1,x2,dataZ["HAsfc9_z"].values)
HAsfc9Model.printParams(x1,x2,dataZ["HAsfc9_z"].values)
The number of levels of the x variables are (2, 11) The standard deviations used for the beta priors are (1.0540496136136044, 2.005676747692769) The standard deviations used for the M12 priors are 0.11328900878057889
pm.model_to_graphviz(HAsfc9Model)
with HAsfc9Model as model:
prior_pred_HAsfc9 = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)
plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_HAsfc9,dataZ["HAsfc9_z"].values,'HAsfc9')
Prior choice is as intended: Broad over the data range.
with HAsfc9Model as model:
trace_HAsfc9 = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.99,random_seed=random_seed)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (10 chains in 10 jobs) NUTS: [HAsfc9_M12_dist, HAsfc9_mu_M12, HAsfc9_b2_dist, HAsfc9_b2_beta, HAsfc9_b1_dist, HAsfc9_b0_dist, HAsfc9_sigmaY, HAsfc9_nuY, HAsfc9_sigma12, HAsfc9_mu_b2, HAsfc9_mu_b1, HAsfc9_mu_b0, HAsfc9_sigma2, HAsfc9_sigma1, HAsfc9_sigma0]
Sampling 10 chains for 1_000 tune and 500 draw iterations (10_000 + 5_000 draws total) took 412 seconds.
with HAsfc9Model as model:
if writeOut:
with open(outPathData + 'model_{}.pkl'.format('HAsfc9'), 'wb') as buff:
pickle.dump({'model': HAsfc9Model, 'trace': trace_HAsfc9}, buff)
with HAsfc9Model as model:
dataTrace_HAsfc9 = az.from_pymc3(trace=trace_HAsfc9)
pm.summary(dataTrace_HAsfc9,hdi_prob=0.95).round(2)
| mean | sd | hdi_2.5% | hdi_97.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| HAsfc9_mu_b0 | -0.07 | 0.80 | -1.65 | 1.43 | 0.01 | 0.01 | 5695.0 | 2359.0 | 5673.0 | 4066.0 | 1.0 |
| HAsfc9_mu_b1[0] | 0.04 | 0.72 | -1.34 | 1.51 | 0.01 | 0.01 | 5050.0 | 2946.0 | 5057.0 | 3783.0 | 1.0 |
| HAsfc9_mu_b1[1] | -0.10 | 0.72 | -1.49 | 1.35 | 0.01 | 0.01 | 4834.0 | 2979.0 | 4834.0 | 3973.0 | 1.0 |
| HAsfc9_mu_b2[0] | -0.26 | 0.66 | -1.56 | 1.03 | 0.01 | 0.01 | 4119.0 | 2839.0 | 4127.0 | 3512.0 | 1.0 |
| HAsfc9_mu_b2[1] | -0.20 | 0.67 | -1.57 | 1.08 | 0.01 | 0.01 | 3738.0 | 3165.0 | 3749.0 | 3996.0 | 1.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| HAsfc9_mu[275] | -0.03 | 0.19 | -0.41 | 0.34 | 0.00 | 0.00 | 4547.0 | 4013.0 | 4560.0 | 4624.0 | 1.0 |
| HAsfc9_mu[276] | 0.06 | 0.24 | -0.39 | 0.52 | 0.00 | 0.00 | 4766.0 | 2887.0 | 4774.0 | 4037.0 | 1.0 |
| HAsfc9_mu[277] | -0.03 | 0.19 | -0.41 | 0.34 | 0.00 | 0.00 | 4547.0 | 4013.0 | 4560.0 | 4624.0 | 1.0 |
| HAsfc9_mu[278] | 0.06 | 0.24 | -0.39 | 0.52 | 0.00 | 0.00 | 4766.0 | 2887.0 | 4774.0 | 4037.0 | 1.0 |
| HAsfc9_mu[279] | -0.03 | 0.19 | -0.41 | 0.34 | 0.00 | 0.00 | 4547.0 | 4013.0 | 4560.0 | 4624.0 | 1.0 |
384 rows × 11 columns
plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_HAsfc9,dataTrace_HAsfc9,'HAsfc9')
/home/bob/.local/lib/python3.8/site-packages/arviz/plots/backends/matplotlib/pairplot.py:212: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid warnings.warn(
with HAsfc9Model as model:
plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_HAsfc9,'HAsfc9')
with HAsfc9Model as model:
plotting_lib.pm.energyplot(trace_HAsfc9)
with HAsfc9Model as model:
posterior_pred_HAsfc9 = pm.sample_posterior_predictive(trace_HAsfc9,samples=numPredSamples,random_seed=random_seed)
/home/bob/.local/lib/python3.8/site-packages/pymc3/sampling.py:1707: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample warnings.warn(
plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_HAsfc9,posterior_pred_HAsfc9,dataZ["HAsfc9_z"].values,'HAsfc9')
plotting_lib.plotLevels(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc9,'HAsfc9',x1,x2)
plotting_lib.plotLevelsStd(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc9,'HAsfc9',x1,x2)
with HAsfc9Model as model:
pm_data_HAsfc9 = az.from_pymc3(trace=trace_HAsfc9,prior=prior_pred_HAsfc9,posterior_predictive=posterior_pred_HAsfc9)
arviz.data.io_pymc3 - WARNING - posterior predictive variable HAsfc9_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.
plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_HAsfc9,'HAsfc9')
plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_HAsfc9,'HAsfc9')
df_hdi_HAsfc9 = plotting_lib.plotTreatmentPosterior(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc9,'HAsfc9',x1,x2)
df_hdi_HAsfc9
| Treatment_i | Treatment_j | hdi_ConfoMap_2.5% | hdi_ConfoMap_97.5% | isSignificant_on_ConfoMap | hdi_Toothfrax_2.5% | hdi_Toothfrax_97.5% | isSignificant_on_Toothfrax | |
|---|---|---|---|---|---|---|---|---|
| 0 | Dry grass | Dry bamboo | 0.000670 | 0.138311 | True | 0.002941 | 0.143668 | True |
| 1 | Dry lucerne | Dry bamboo | -0.035688 | 0.084977 | False | -0.042002 | 0.076006 | False |
| 2 | Dry lucerne | Dry grass | -0.117764 | 0.031872 | False | -0.125602 | 0.021941 | False |
| 3 | BrushNoDirt | BrushDirt | -0.231361 | 0.146793 | False | -0.063532 | 0.221075 | False |
| 4 | Control | BrushDirt | -0.159309 | 0.177833 | False | -0.020488 | 0.239702 | False |
| 5 | Control | BrushNoDirt | -0.118868 | 0.215704 | False | -0.121261 | 0.192874 | False |
| 6 | RubDirt | BrushDirt | -0.223620 | 0.119671 | False | -0.093554 | 0.146212 | False |
| 7 | RubDirt | BrushNoDirt | -0.184998 | 0.151141 | False | -0.204337 | 0.093804 | False |
| 8 | RubDirt | Control | -0.211662 | 0.080935 | False | -0.212686 | 0.056618 | False |
| 9 | Clover+dust | Clover | -0.147903 | 0.136656 | False | -0.081257 | 0.155239 | False |
| 10 | Grass | Clover | -0.251859 | 0.071005 | False | -0.126162 | 0.116917 | False |
| 11 | Grass | Clover+dust | -0.231351 | 0.051562 | False | -0.160442 | 0.079536 | False |
| 12 | Grass+dust | Clover | -0.310090 | 0.004148 | False | -0.202760 | 0.043221 | False |
| 13 | Grass+dust | Clover+dust | -0.272959 | -0.010159 | True | -0.240244 | -0.005423 | True |
| 14 | Grass+dust | Grass | -0.211711 | 0.091591 | False | -0.206979 | 0.039560 | False |
if writeOut:
df_hdi_HAsfc9.to_csv(outPathData+ 'hdi_{}.csv'.format('HAsfc9'))
with pm.Model() as model:
HAsfc81Model = TwoFactorModel('HAsfc81',x1,x2,dataZ["HAsfc81_z"].values)
HAsfc81Model.printParams(x1,x2,dataZ["HAsfc81_z"].values)
The number of levels of the x variables are (2, 11) The standard deviations used for the beta priors are (1.0444803217312628, 1.586983908089902) The standard deviations used for the M12 priors are 0.09499285587637817
pm.model_to_graphviz(HAsfc81Model)
with HAsfc81Model as model:
prior_pred_HAsfc81 = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)
plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_HAsfc81,dataZ["HAsfc81_z"].values,'HAsfc81')
Prior choice is as intended: Broad over the data range.
with HAsfc81Model as model:
trace_HAsfc81 = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.99,random_seed=random_seed)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (10 chains in 10 jobs) NUTS: [HAsfc81_M12_dist, HAsfc81_mu_M12, HAsfc81_b2_dist, HAsfc81_b2_beta, HAsfc81_b1_dist, HAsfc81_b0_dist, HAsfc81_sigmaY, HAsfc81_nuY, HAsfc81_sigma12, HAsfc81_mu_b2, HAsfc81_mu_b1, HAsfc81_mu_b0, HAsfc81_sigma2, HAsfc81_sigma1, HAsfc81_sigma0]
Sampling 10 chains for 1_000 tune and 500 draw iterations (10_000 + 5_000 draws total) took 718 seconds.
with HAsfc81Model as model:
if writeOut:
with open(outPathData + 'model_{}.pkl'.format('HAsfc81'), 'wb') as buff:
pickle.dump({'model': HAsfc81Model, 'trace': trace_HAsfc81}, buff)
with HAsfc81Model as model:
dataTrace_HAsfc81 = az.from_pymc3(trace=trace_HAsfc81)
pm.summary(dataTrace_HAsfc81,hdi_prob=0.95).round(2)
| mean | sd | hdi_2.5% | hdi_97.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| HAsfc81_mu_b0 | -0.06 | 0.81 | -1.66 | 1.48 | 0.01 | 0.01 | 5438.0 | 3055.0 | 5456.0 | 4451.0 | 1.0 |
| HAsfc81_mu_b1[0] | 0.01 | 0.72 | -1.40 | 1.39 | 0.01 | 0.01 | 4786.0 | 2800.0 | 4784.0 | 3776.0 | 1.0 |
| HAsfc81_mu_b1[1] | -0.08 | 0.74 | -1.54 | 1.35 | 0.01 | 0.01 | 4760.0 | 2996.0 | 4760.0 | 3738.0 | 1.0 |
| HAsfc81_mu_b2[0] | -0.36 | 0.67 | -1.69 | 0.90 | 0.01 | 0.01 | 3853.0 | 3574.0 | 3866.0 | 3483.0 | 1.0 |
| HAsfc81_mu_b2[1] | -0.40 | 0.67 | -1.68 | 0.92 | 0.01 | 0.01 | 4465.0 | 3899.0 | 4474.0 | 3963.0 | 1.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| HAsfc81_mu[275] | -0.16 | 0.17 | -0.45 | 0.16 | 0.00 | 0.00 | 5195.0 | 5195.0 | 5325.0 | 4756.0 | 1.0 |
| HAsfc81_mu[276] | -0.10 | 0.21 | -0.44 | 0.34 | 0.00 | 0.00 | 4753.0 | 4660.0 | 5116.0 | 4932.0 | 1.0 |
| HAsfc81_mu[277] | -0.16 | 0.17 | -0.45 | 0.16 | 0.00 | 0.00 | 5195.0 | 5195.0 | 5325.0 | 4756.0 | 1.0 |
| HAsfc81_mu[278] | -0.10 | 0.21 | -0.44 | 0.34 | 0.00 | 0.00 | 4753.0 | 4660.0 | 5116.0 | 4932.0 | 1.0 |
| HAsfc81_mu[279] | -0.16 | 0.17 | -0.45 | 0.16 | 0.00 | 0.00 | 5195.0 | 5195.0 | 5325.0 | 4756.0 | 1.0 |
384 rows × 11 columns
plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_HAsfc81,dataTrace_HAsfc81,'HAsfc81')
/home/bob/.local/lib/python3.8/site-packages/arviz/plots/backends/matplotlib/pairplot.py:212: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid warnings.warn(
with HAsfc81Model as model:
plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_HAsfc81,'HAsfc81')
with HAsfc81Model as model:
plotting_lib.pm.energyplot(trace_HAsfc81)
with HAsfc81Model as model:
posterior_pred_HAsfc81 = pm.sample_posterior_predictive(trace_HAsfc81,samples=numPredSamples,random_seed=random_seed)
/home/bob/.local/lib/python3.8/site-packages/pymc3/sampling.py:1707: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample warnings.warn(
plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_HAsfc81,posterior_pred_HAsfc81,dataZ["HAsfc81_z"].values,'HAsfc81')
plotting_lib.plotLevels(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc81,'HAsfc81',x1,x2)
with HAsfc81Model as model:
pm_data_HAsfc81 = az.from_pymc3(trace=trace_HAsfc81,prior=prior_pred_HAsfc81,posterior_predictive=posterior_pred_HAsfc81)
arviz.data.io_pymc3 - WARNING - posterior predictive variable HAsfc81_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.
plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_HAsfc81,'HAsfc81')
plotting_lib.plotLevelsStd(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc81,'HAsfc81',x1,x2)
plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_HAsfc81,'HAsfc81')
df_hdi_HAsfc81 = plotting_lib.plotTreatmentPosterior(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc81,'HAsfc81',x1,x2)
df_hdi_HAsfc81
| Treatment_i | Treatment_j | hdi_ConfoMap_2.5% | hdi_ConfoMap_97.5% | isSignificant_on_ConfoMap | hdi_Toothfrax_2.5% | hdi_Toothfrax_97.5% | isSignificant_on_Toothfrax | |
|---|---|---|---|---|---|---|---|---|
| 0 | Dry grass | Dry bamboo | -0.009401 | 0.101134 | False | 0.018900 | 0.136196 | True |
| 1 | Dry lucerne | Dry bamboo | -0.030370 | 0.087266 | False | -0.001664 | 0.109981 | False |
| 2 | Dry lucerne | Dry grass | -0.073950 | 0.050741 | False | -0.087209 | 0.039690 | False |
| 3 | BrushNoDirt | BrushDirt | -0.518195 | 0.206099 | False | -0.147466 | 0.286371 | False |
| 4 | Control | BrushDirt | -0.376578 | 0.303712 | False | -0.030208 | 0.268768 | False |
| 5 | Control | BrushNoDirt | -0.097626 | 0.347835 | False | -0.157188 | 0.295603 | False |
| 6 | RubDirt | BrushDirt | -0.493950 | 0.236296 | False | -0.128138 | 0.198331 | False |
| 7 | RubDirt | BrushNoDirt | -0.231775 | 0.307111 | False | -0.273242 | 0.218386 | False |
| 8 | RubDirt | Control | -0.321828 | 0.132966 | False | -0.274491 | 0.088330 | False |
| 9 | Clover+dust | Clover | -0.455539 | 0.323286 | False | -0.516906 | 0.156140 | False |
| 10 | Grass | Clover | -0.656279 | 0.348255 | False | -0.691228 | 0.029980 | False |
| 11 | Grass | Clover+dust | -0.454243 | 0.139046 | False | -0.317566 | 0.103036 | False |
| 12 | Grass+dust | Clover | -1.007356 | -0.112084 | True | -0.965316 | -0.269417 | True |
| 13 | Grass+dust | Clover+dust | -0.689484 | -0.292849 | True | -0.567226 | -0.229459 | True |
| 14 | Grass+dust | Grass | -0.708570 | -0.038109 | True | -0.516849 | -0.061466 | True |
if writeOut:
df_hdi_HAsfc81.to_csv(outPathData+ 'hdi_{}.csv'.format('HAsfc81'))
For e.g. the pair Clover+Dust vs. Clover shows an unexpected bimodal distribution of the contrast.
We now examine the traces carefully to exclude errors in the sampling:
Get the traces on ConfoMap of the interactions
m12_clover = trace_HAsfc81['HAsfc81_M12'][:,0,2]
m12_cloverDust = trace_HAsfc81['HAsfc81_M12'][:,0,3]
diff_m12 = m12_cloverDust-m12_clover
Get the traces of the treatments:
b2_clover = trace_HAsfc81['HAsfc81_b2'][:,2]
b2_cloverDust = trace_HAsfc81['HAsfc81_b2'][:,3]
diff_b2 = b2_cloverDust-b2_clover
Look at all the pairs
sns.pairplot(data=pd.DataFrame.from_dict(
{'m12_clover':m12_clover,'m12_cloverDust':m12_cloverDust,'diff_m12':diff_m12,
'b2_clover':b2_clover,'b2_cloverDust':b2_cloverDust,'diff_b2':diff_b2
}
));
Let's have that of the differences again
sns.jointplot(x=diff_b2,y=diff_m12);
We see two sets that are distributed along parallel lines in the scatter plot.
This means that the model estimates two subsets of possible differences.
However, when looking at the raw data at 'analysis/plots/SSFA_Sheeps_plot.pdf' one can see that the distribution of values for HAsfc81 on Clover (Sheep) measured by ConfoMap appear to have a bimodal distribution.
Thus, in combination with the chosen uninformative priors, the model correctly describes the distribution as bimodal.
In summary, we see no issue with the modeling and sampling.
Set the surface parameters for every treatment dataframe:
df_hdi_Asfc["SurfaceParameter"] = "Asfc"
df_hdi_HAsfc9["SurfaceParameter"] = "HAsfc9"
df_hdi_HAsfc81["SurfaceParameter"] = "HAsfc81"
df_hdi_R["SurfaceParameter"] = "R²"
df_hdi_epLsar["SurfaceParameter"] = "epLsar"
df_hdi_total = pd.concat([df_hdi_epLsar,df_hdi_R,df_hdi_Asfc,df_hdi_HAsfc9,df_hdi_HAsfc81],ignore_index=True)
Show the treatment pairs and surface parameters where the softwares differ
df_summary = df_hdi_total[df_hdi_total.isSignificant_on_ConfoMap != df_hdi_total.isSignificant_on_Toothfrax][["Treatment_i","Treatment_j","SurfaceParameter","isSignificant_on_ConfoMap","isSignificant_on_Toothfrax","hdi_ConfoMap_2.5%","hdi_ConfoMap_97.5%","hdi_Toothfrax_2.5%","hdi_Toothfrax_97.5%"]]
df_summary
| Treatment_i | Treatment_j | SurfaceParameter | isSignificant_on_ConfoMap | isSignificant_on_Toothfrax | hdi_ConfoMap_2.5% | hdi_ConfoMap_97.5% | hdi_Toothfrax_2.5% | hdi_Toothfrax_97.5% | |
|---|---|---|---|---|---|---|---|---|---|
| 2 | Dry lucerne | Dry grass | epLsar | True | False | 0.000228 | 0.001524 | -0.000167 | 0.001116 |
| 10 | Grass | Clover | epLsar | False | True | -0.000153 | 0.002973 | 0.000293 | 0.003795 |
| 18 | BrushNoDirt | BrushDirt | R² | True | False | -0.002096 | -0.000672 | -0.000420 | 0.000635 |
| 22 | RubDirt | BrushNoDirt | R² | True | False | 0.000181 | 0.002386 | -0.000556 | 0.000619 |
| 24 | Clover+dust | Clover | R² | True | False | -0.001005 | -0.000290 | -0.000138 | 0.000401 |
| 26 | Grass | Clover+dust | R² | False | True | -0.000639 | 0.000580 | -0.000792 | -0.000199 |
| 28 | Grass+dust | Clover+dust | R² | False | True | -0.001314 | 0.000824 | -0.000515 | -0.000045 |
| 33 | BrushNoDirt | BrushDirt | Asfc | True | False | 0.013059 | 7.915624 | -2.631889 | 6.887819 |
| 35 | Control | BrushNoDirt | Asfc | True | False | 0.022112 | 14.415237 | -0.616464 | 13.182578 |
| 36 | RubDirt | BrushDirt | Asfc | False | True | -0.259987 | 20.831039 | 1.502426 | 13.522475 |
| 60 | Dry grass | Dry bamboo | HAsfc81 | False | True | -0.009401 | 0.101134 | 0.018900 | 0.136196 |
if writeOut:
df_summary.to_csv(outPathData+ 'summary.csv')
!jupyter nbconvert --to html Statistical_Model_TwoFactor.ipynb
[NbConvertApp] Converting notebook Statistical_Model_TwoFactor.ipynb to html [NbConvertApp] Writing 19659153 bytes to Statistical_Model_TwoFactor.html
!jupyter nbconvert --to markdown Statistical_Model_TwoFactor.ipynb
[NbConvertApp] Converting notebook Statistical_Model_TwoFactor.ipynb to markdown [NbConvertApp] Support files will be in Statistical_Model_TwoFactor_files/ [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Making directory Statistical_Model_TwoFactor_files [NbConvertApp] Writing 94860 bytes to Statistical_Model_TwoFactor.md